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We compute by means of Quantum Monte Carlo simulations the equation of state of bulk solid 
parahydrogen extrapolated to zero temperature, up to a pressure of ~ 2 MBar. We compare the 
equation of state yielded by three different pair potentials, namely the Silvera-Goldman, Buck and 
one recently proposed by Moraldi, modified at short distances to include a repulsive core, missing 
in the originally proposed potential. The Moraldi pair potential yields an equation of state in very 
good agreement with experiment at megabar pressures, owing to its softer core, and is at least as 
accurate as the SG or the Buck at saturated vapour pressure. Estimates for the experimentally 
measurable kinetic energy per molecule are provided for all pair potentials. 



I. INTRODUCTION 

Hydrogen is the simplest and most abundant element 
in the universe. Achieving predictive knowledge of its 
equation of state in a broad range of thermodynamic 
conditions remains a worthwhile theoretical goal, due to 
its relevance to a wide variety of physical systems, of- 
ten with possible technological applications,'^^ and also 
because it provides a cogent test of the most advanced 
computational many-body techniques. 

The equation of state (EOS) of hydrogen, in essen- 
tially all phases of interest, is a fully quantum- mechanical 
many-body problem. Theoretical calculations broadly 
fall into two different categories; in ah initio studies, 
the mathematical model explicitly takes into account 
the elementary constituents (namely electrons and pro- 
tons), which interact via the electrostatic Coulomb po- 
tential. Calculations are carried out within the frame- 
work of Density Functional Theory (DFTpHSl 

or Quan- 
tum Monte Carlo (QMC).li2Ml Qn the other hand, in 
molecular condensed phases one often adopts the Born- 
Oppenheimer approximation, allowing one to regard in- 
dividual hydrogen molecules as elementary particles, and 
describe their interaction by means of a static potential, 
the simplest choice being that of a central pair potential. 

An approach based on a pair potential only depen- 
dent on the distance between two molecules, obviously 
cannot describe any process involving electronic trans- 
fer, nor the energy contribution of interactions involv- 
ing, e.g., triplets of molecules, nor any effect arising from 
the non-spherical character of the inter-molecular inter- 
action; the importance of all of these physical mecha- 
nisms generally increases with the thermodynamic pres- 
sure of the condensed phase under investigation. How- 
ever, there are distinct advantages to the use of static 
pair-wise potentials, mainly that it is typically computa- 
tionally much faster and conceptually simpler than ah 
initio approaches. Furthermore, if pair potentials are 
used, thermodynamic properties of molecular hydrogen 
in the condensed phase can typically be computed es- 
sentially without any uncontrolled approximation, e.g., 
by means of QMC simulations. Thus, it is desirable to 
develop pair potentials affording a reasonably accurate, 
quantitative description of the condensed phase of hydro- 
gen in broad ranges of thermodynamic conditions. 

The Silvera-Goldman model pair potentiaP^ is ar- 



guabl y t he most commonly adopted, and has been 
showrpSfo afford a quantitative description of the the low 
temperature, equilibrium solid phase of parahydrogen 
(the only species discussed here). Another popular pair 
potential is the Buck,'i2H2i] which is very similar to the 
Silvera-Goldman, featuring a slightly deeper attractive 
well. Both potentials yield an EOS in reasonable agree- 
ment with experiment at moderate pressure (< 25 kPa). 
At higher pressure, however, they become increasingly in- 
accurate, leading to an overstimation of the pressure due 
to the excessive "stiffness" of their repulsive core at short 
inter- molecular separation (below ~ 2.8 A). This problem 
is also present in theoretical calculations of the EOS of 
condensed helium at high pressure, based on pair poten- 
tials; in that context, it is known that better agreement 
with the experimental EOS can be obained on including 
three-body terms, whose overall effect is that of softening 
the repulsive core of the pair-wise interaction.'^^ltMi Still, 
given the significant computational overhead involved in 
the inclusion of three-body terms, the question remains 
of whether a modified effective pair potential could offer 
more satisfactory agreement between theory and experi- 
ment at high (e.g., megabar) pressure. 

To this aim, Moraldi recently proposecP^ a modified 
effective pair potential, with a softened repulsive core, 
designed to afford greater agreement between the theo- 
retically predicted values of the pressure as a function of 
the density in the limit of low temperature (T — ^ 0). Re- 
sults of perturbative calculations carried out in Ref. [55] 
show remarkable agreement with experimental data of 
the pressure computed with the modified pair potential, 
in an extended range of density (up to 0.24 A~^, which 
corresponds to a pressure in the megabar range). 

In this work, the T = EOS of solid parahydrogen is 
computed by means of first principle QMC simulations, 
in a wider density range with respect to Ref. 1251 namely 
up to a density of 0.273 A~'^, for a model of condensed 
parahydrogen based on three different pair potentials, 
namely a slightly modified version of the Moraldi po- 
tential (see below), the Silvera-Goldmann and the Buck. 
Our purpose is that of providing a numerically accurate, 
non-perturbative check of the pressure computed in Ref. 
1251 and to offer a cogent comparison of the performance 
of the three potentials. We also compute the total and 
kinetic energy per molecule; the latter is experimentally 
measurable by means of neutron scattering, and is espe- 
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cially sensitive to the detailed features of the repulsive 
core of the interactionl^^l 

We find that the recently proposed pair potential in- 
deed yields pressure estimates in much better agreement 
with experiment than the Silvera-Goldman and Buck po- 
tential, especially above 0.14 A^'^. In fact, the agree- 
ment between the values of the pressure computed with 
the modified Moraldi potential and the experimentally 
measured ones remains excellent up to ~ 180 GPa, the 
highest pressure for which experimental data are avail- 
able for solid H2. At the equilibrium density, i.e., zero 
pressure at r=0, the Moraldi potential yields energy esti- 
mates comparable to those of the Silvera-Goldman. Ki- 
netic energy estimates at high density computed with 
the Moraldi potential are in significant disaccord with 
(considerably lower than) those furnished by the Silvera- 
Goldman and Buck potentials, a fact ascribable to the 
softer repulsive core of the Moraldi pair interaction. A 
measurement of the kinetic energy will therefore furnish 
additional important insight on whether the Moraldi po- 
tential not only yields an EOS in quantitatively close 
agreement with experiment, but also affords a generally 
more accurate, quantitatively physical description of the 
system than the Silvera-Goldman or Buck model inter- 
actions. 

In the next section, the microscopic model underlying 
the calculation, specifically the pair potentials utilized 
are illustrated; we briefly outline the basic methodology 
utilized and offer relevant computational details in Sec. 
|III| In Sec. |IV| a thorough illustration of the results ob- 
tained in this work is provided; finally, a general discus- 
sion is offered, and conclusions outlined, in Sec. [V] 



principle, preventing the electronic clouds of different 
molecules from spatially overlapping, whereas Vatt(r) 
represents the long-range Van der Waals attraction of 
mutually induced molecular electric dipoles. The func- 
tion fc {r) has the purpose of removing the divergence of 
Vattir) as r — > 0. The detailed forms of these functions, 
as well as the values of the parameters for the two po- 
tentials, are furnished in Refs. fTTIandfTOl 

The Moraldi potential, as proposed in Ref. is ex- 
pressed as follows: 



VMir) = Vrepir).frep{r) - Vattir) fc (r) 



(3) 



where the functions Vrep, Vatt and fc are those of the 
Silvera-Goldman potential, whereas frep{f) is a function 
of the same form of fc, softening the repulsion at short 
distances contained in Vrep- As shown in Ref. 1251 the 
potential Vm('') displays a much slower growth as r — > 
than the Silvera-Goldman. In fact, the function VMii") 
stops growing altogether at r ^ r-o = 1.28 A, and actu- 
ally decreases all the way to zero as r — >■ (see Figure [I]). 
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II. MODEL 

Our system of interest is modeled as an ensemble of N 
parahydrogen molecules, regarded as point particles of 
spin zero, enclosed in a vessel of volume i7, shaped as a 
parallelepiped, with periodic boundary conditions in all 
directions. The sides of the parallelepiped are chosen to 
fit a crystalline sample of solid parahydrogen (hexagonal 
close-packed structure) . The quantum- mechanical many- 
body Hamiltonian is the following: 

JV 

H = -\Y,vUY.^{ni) (1) 

1=1 i<j 

Here, A=12.031 kA^, while V is the potential describing 
the interaction between two molecules, only depending 
on their relative distance. Because this work is specif- 
ically about comparing different potentials, we provide 
here some basic details of all potentials used. The Silvera- 
Goldman {Vsg) and Buck {Vb) potentials have the form 

V{r) = Vrep{r)-Vatt{r)fc{r) (2) 

where Vrepir) describes the repulsion of two molecules 
at short distances, arising mainly from Pauli exclusion 
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FIG. 1: Color online. The Moraldi pair potential VA,f(r) as 
proposed in Ref. [55] (solid line) and the modified one U (r) 
used in the present work (dashed line). 

This unphysical feature of Vjvf (r) is of little consequences 
in calculations at low density, but must be corrected in 
order to prevent particles from "tunnelling" across the 
potential barrier and falling on top of one another, some- 
thing observed in computer simulations making use of 
Vm {t) at densities for which the inter- molecular distance 
is < 2 A. Thus, in this work we made use of the modified 
potential C/(r), defined as follows: 

U{r) = VM{r)Q{r - d) + Vrepir)e{d - r) (4) 

where d—A a.u. (i.e., 2.1167 A), 0(a;) is the Heaviside 
function and Vrep{r) has the same form of Vrep, but the 
two parameters (usually referred to as a and /3) upon 
which it depends are adjusted to match Vuir) and its 
first derivative at r — d. This modified potential, utilized 
in our simulations, is shown alongside VM{r) in Figure [T] 
It does not have the unphysical feature of VM{r) at short 
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distances discussed above, indeed it increases monotoni- 
cally as r — )■ 0. 

The above modification of tlie Moraldi potential is ob- 
viously one of the many that are possible; it has the ad- 
vantage of being easy to implement. The detailed form of 
the repulsive core of the potential at molecular separation 
shorter than ~ 1.5 A turns out not to affect the calcula- 
tion significantly, in that particles seldom find themselves 
at such distances, in the range of pressure explored here. 
In other words, a repulsive core is needed, in order to 
avoid the unphysical behavior described above, but its 
precise form is not crucial, in the range of density con- 
sidered in this work. 



by illustrating results of QMC simulations at the satura- 
tion density po = 0.0261 A~^. Specifically, we compare 
the pressure and energy estimates yielded by the Moraldi 
potential (modified as explained above), which was not 
purposefully designed to describe the low-density equi- 
librium phase, to those yielded by the Silvera-Goldman 
and Buck potentials, as well as to experimental data. 
The purpose of such a comparison is that of ascertain- 
ing whether the softening of the Silvera-Goldman repul- 
sive core, which is at the basis of the Moraldi potential, 
might worsen the agreement with experiment at low den- 
sity, while possibly improving it in the megabar pressure 
range. 



III. CALCULATION 

The T = (ground state) equation of state of solid 
parahydrogen, modelled by the many-body Hamiltonian 
([T]) with the three different potentials Vsg, Vb and U 
described above, was computed by means of numerical 
QMC simulations based on the continuous-space Worm 
Algorithm (WA).'23l2ll Specifically, we estimate the pres- 
sure P as a function of the density p, in the interval 
0.0763 < P < 0.273 A"^, for a finite system com- 
prising A'^=216 molecules. In order to provide ground 
state estimates, we perform simulations at different low 
temperatures and extrapolate the results to T = 0. We 
generally find that results obtained at T = 4 K are in- 
distinguishable from the extrapolated ones, within our 
quoted statistical uncertainties. We use the standard 
virial estimator^ for the pressure, and estimate the con- 
tribution from particles lying outside the main simula- 
tion cell by radially integrating the quantity r{dV/dr), 
approximating the pair correlation function g{r) to 1 
(we do the same for the potential energy). ■^'^ The QMC 
results presented here are obtained using a time step 
T = 5 X 10"'' K"'^, as we observed that the estimates 
yielded by such a choice of time step coincide with those 
extrapolated to r = 0, within statistical uncertainties. 

It should be mentioned that in principle parahydrogen 
molecules must be regarded as indistinguishable parti- 
cles of spin zero, and thus obeying Bosc statistics. The 
WA explicitly allows for permutations of identical par- 
ticles, and quantum-mechanical exchanges can be im- 
portant in parahydrogen, in specific circumstances. For 
example, they are predicted to underlie superfiuidity 
at low te mpera ture in small clusters (less than thirty 
molecules) pTT32| g^j-^^ their effect is measurable in the mo- 
mentum distribution of the liquid near melting.'^ In the 
solid phase, however, in the range of density explored 
here, exchanges are practically absent, i.e., particles can 
be regarded for practical purposes as distinguishable. 



IV. RESULTS 

Although the focus of our work is directed to the 
ground state EOS of solid parahydrogen at relatively high 
(megabar) pressure, we begin the discussion of the results 





This work 


Ref. Lia 


Moraldi 


-88.5(2) 




Silvera-Goldman 


-88.1(1) 


-87.90(2) 


Buck 


-93.8(1) 


-93.87(2) 



TABLE I: Ground state energy per molecule (in K) for solid 
parahydrogen in the hep phase at the saturation density po = 
0.0261 A"'^, computed in this work and in Ref. [l8]for various 
intermolecular potentials. Statistical errors, in parentheses, 
are on the last digit. 

Table |l] shows numerical estimates of the ground state 
energy per molecule for solid parahydrogen at saturation 
density. The estimates obtained here, obtained by ex- 
trapolating all the way down to T = results obtained 
in the range 1 K < T < 4 K, are in agreement with 
those obtained in Ref. [15] using ground state Diffusion 
Monte Carlo (DMC) simulations. Experimentally, the 
ground state energy per molecule has been inferred in 
different ways, and the agreement between the various 
estimatea^SEZl is not perfect; quoted values range from 
-89.9 K of Ref. f^to -93.5 K of Refs. and [371, it is 
clear from the data reported in the table that the Moraldi 
potential affords an accuracy comparable to that of the 
other two potentials. 

We compare in figure [2] our computed values of the 
pressure for hep parahydrogen in the T limit, in a 
density range up to 0.273 A~'^, which corresponds to a 
pressure of about 1.8 Mbars. We also include for com- 
parison a fit to the most recenlPsl compilation of experi- 
mental (x-ray) measurement a'^'^^ ' "^"' of the EOS obtained 
with the Vinet equation of state,!^ 



3/fo(l-rf) 



(K-m-d) 



(5) 



where d = {pa/ pY^^, with parameters Kq ~ 0.20(1) GPa 
and K'^ = 6.84(7). The computed values of the pressure 
are also reported in Table [lT| 

The agreement between the experimental values of 
the pressure, and those obtained by simulation using 
the Moraldi potential is excellent, considering the rela- 
tively simple form of the pair potential itself. At the 
highest compression considered here, namely 0.273 A~'^, 
the Silvera-Goldman potential overestimate the pressure 
by over a factor two. The Buck potential does only 



4 



ci3 
O 



3 

1/3 



400 ■ 



300 ■ 



200 ■ 



100 ■ 



• Buck 

♦ SG 

■ Moraldi 
-- Experiment 



-■ < 



0.08 0.12 0.16 0.20 0.24 0.28 
Density (A^) 



FIG. 2: Color online. Equation of state of solid parahy- 
drogen at T=0, computed by Quantum Monte Carlo simu- 
lations using the Silvera- Goldman (diamonds), Buck (circles) 
and Moraldi (boxes) pair potentials. Dashed line represents 
a fit to experimental data from Ref. [SSj obtained with the 
Vinet equation of state (see text). Combined statistical and 
systematic errors affecting our computed pressures are at the 
most of the order of 0.1% of the pressure. 



Density Moraldi Buck 



SG 



0.076 
0.098 
0.129 
0.174 
0.243 
0.273 



5.11 X 10* 6.09 X 10" 5.92 x 10* 

1.17 X 10^ 1.54 X 10^ 1.52 x 10^ 

2.61 X 10^ 3.87 X 10^ 3.93 x 10"^ 

5.93 X 10^ 9.74 x 10^ 1.04 x 10*^ 

1.34 X 10® 2.47 X lO'' 2.79 x lO'' 

1.80 X 10** 3.40 X 10® 3.94 x 10® 



TABLE II: Pressure (in bars) calculated for different densities 
(in A""*) in the T — > limit, using Buck, Silvera-Goldman 
and Moraldi potentials. Combined statistical and systematic 
errors affecting our QMC simulations are estimated to be of 
the order of 0.1% or less of the quoted values. 



marginally better than the Silvera-Goldman. The lower 
values of the pressure yielded by the Moraldi potential are 
a direct consequence of its softer repulsive core. It should 
be noted that the values in Table |lT] are in quantitative 
agreement with those obtained by Moraldi in Ref. 
based on a perturbative calculation. Perhaps even more 
significant is the fact that agreement remains excellent 
even above 110 GPa, i.e., a compression at which solid 
hydrogen undergoes an orientational phase transition.^ 
One would expect that in the orientationally ordered 
phase, the approximation of a spherically symmetric po- 
tential would be particularly problematic. 

The fact that the Moraldi potential is not as stiff as 
the Sivera-Goldman or Buck potentials at short distances 
is also reflected on the predicted structural properties of 
the crystal. In particular, the local environment experi- 
enced by a single atom is captured by the static structure 
factor, which of course can be probed by X-ray diffrac- 
tion. Its Fourier transform, the pair correlation function 
g{r), is easily accessible by simulation. Results for two 




FIG. 3: Color online. Pair correlation function g{r) for hep 
parahydrogen at T—4: K, computed by simulation at the two 
densities p=0.0763 A"^ and p=0.2428 A"^, with the Moraldi 
and Buck potentials. The results for the Silvera-Goldman 
potential are similar to those for the Buck. 



different densities are shown in Fig. [Sj where the pair 
correlation functions computed for the Buck and Moraldi 
potentials are compared. It should be noted that the SG 
and Buck potentials yield similar results for this quan- 
tity (the main peak of the g{r) is ^ 10% higher with the 
SG potential, at the higher density). While at the lower 
density, which corresponds to a pressure around 5 kilo- 
bars, the pair correlation functions yielded by the two 
potentials are comparable, at the higher density (which 
corresponds to a pressure close to one Megabar), there 
is a clear difference between the two results. Specifically, 
the pair correlation function computed with the Buck po- 
tential displays considerably higher peaks and generally 
much sharper features and an overall more classical struc- 
ture, compared to that obtained with the softer Moraldi 
potential, which is smoother, and displays a markedly 
more quantum-mechanical character. 

While the improvement afforded by the Moraldi po- 
tential on the calculation of the EOS is certainly remark- 
able, especially considering the relatively wide range of 
pressure for which it is observed, some care should be 
exercised if a comprehensive assessment is sought of the 
microscopic physical description yielded by a given pair 
potential. It is possible to obtain, by means of an ef- 
fective pair potential, an EOS in close agreement with 
experiment, but at the cost of worsening the agreement 
for other quantities. 

In particular, the kinetic energy per particle is particu- 
larly sensitive to the detailed shape of the repulsive core 
of the intermolecular potential at short distances (see, for 
instance, discussion in Ref. [22]) . and it can also be mea- 
sured experimentally by neutron scattering, as the second 
moment of the single-particle momentum distribution.!^ 
Thus, a cogent test of the relative accuracy of the new po- 
tential with respect to the existing ones could be offered 
by measurements of the kinetic energy per particle, for 
which we have obtained ground state estimates as well, 
for all three potentials considered here. The results are 
displayed in Fig. |4]and detailed in Table [TTTl 

At the saturation density po all potentials yield essen- 
tially the same kinetic energy per molecule, whereas at 
higher density the softer core of the Moraldi potential 
results in a considerably lower kinetic energy, especially 
at the highest compression. It should be noted, in any 
case, that the kinetic energy contribution to the pressure 
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FIG. 4: Color online. Kinetic energies per molecule at r=0, 
computed by simulations at different densities for the three 
potentials considered here 



Density Moraldi SG Buck RefESl 



0.0261 


69 


70 


71 




0.0763 


341 


406 


411 


405 


0.0980 


440 


576 


570 


531 


0.1288 


558 


817 


786 


678 


0.1739 


701 


1170 


1060 


834 


0.2428 


872 


1650 


1390 


914 


0.2730 


927 


1855 


1523 





TABLE III: Kinetic energy per molecule (in K) in the T — >■ 
limit, at various densities (in A~"^), obtained using each 
of the three pair potentials considered in this work in the 
QMC simulation. Combined statistical and systematic errors 
affecting our results are estimated at one percent or less of 
the quoted values. Also reported are the estimates from Ref. 
1251 in the rightmost column. 



in this work with those reported in Ref. shows that 
the perturbative approach leads to a slight overestima- 
tion, but otherwise furnishes reasonably accurate results 
for the problem of interest. 



V. CONCLUSIONS 

In this work we have carried out a numerical analysis 
of a new intermolecular potential proposed by Moraldi, 
aimed at reproducing the experimental equation of 
state of solid parahydrgen at low temperature, up to a 
pressure of 1.8 Mbars. The potential considered here is 
a modification of the one proposed in Ref. [211 to which 
a repulsive core has been added at short distances, to 
prevent unphysical behavior at the highest densities 
considered. This potential has a considerably softer core 
than the Silvera-Goldman and Buck ptential, resulting 
in enhanced quantum effects and lower values of the 
pressure. We have carried out Quantum Monte Carlo 
simulation of hep parahydrogen in the T -> limit, and 
compared estimates yielded by the new potential, as 
well as the Silvera-Goldman and Buck potentials, for 
the pressure and for the kinetic energy per molecule, 
as a function of density. We have also compared the 
numerical results to the most recent experimental 
measurements of the equation of state, and found that 
the new potential yields excellent agreement with ex- 
periment, agreement which extends all the way down to 
low density (saturation). We also furnish estimates for 
the kinetic energy per molecule, whose measurement by 
neutron diffraction will provide another independent, im- 
portant assessment of the reliability of the new potential. 
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